library(readr)
setwd('F:/Rproject/baseline')
mydata <- read.csv("hypoglycemia_end.csv")
# mydata <- post_thrombotic_syndrome
dim(mydata)

#删除有缺失值的行
mydata <- na.omit(mydata)
dim(mydata)
View(mydata)
names(mydata)
str(mydata)
head(mydata)
summary(mydata)

#将分类变量转换为因子
mydata$marital_status <- factor(mydata$marital_status,
                                levels = c(1,2,3,4),
                                labels = c('married','unmarried','divorce','widow'))
str(mydata)
table(mydata$marital_status)

#拆分数据
dev <- mydata[mydata$dataset==1,]
vad <- mydata[mydata$dataset==0,]

#循环创建单因素回归模型
uni_glm_model <- function(x){
  #paste0()类似于字符串拼接
  #hypoglycemia==1表示低血糖
  #as.formula()将字符串转换为公式
  FML <- as.formula(paste0("hypoglycemia==1~",x))
  glm1 <- glm(FML,data = dev,family = binomial)
  glm2 <- summary(glm1)
  
  #计算获取需要的指标
  OR <- round(exp(coef(glm1)),2)
  SE <- round(glm2$coefficients[,2],3)
  CI2.5 <- round(exp(coef(glm1)-1.96*SE),2)
  CI97.5 <- round(exp(coef(glm1)+1.96*SE),2)
  CI<-paste0(CI2.5,'-',CI97.5)
  B<-round(glm2$coefficients[,1],3)
  Z<-round(glm2$coefficients[,3],3)
  P<-round(glm2$coefficients[,4],3)
  
  
  
  #将计算出来的指标制作为数据框
  uni_glm_model<-data.frame('characteristics'=x,
                            'B'=B,
                            'SE'=SE,
                            'OR'=OR,
                            'CI'=CI,
                            'Z' =Z,
                            'P'=P)[-1,]
  
  return(uni_glm_model)
}
#指定参与循环分析的若干自变量X
variable.names <- colnames(dev[c(4:40)])
variable.names

#运行上面自定义批量执行函数
#lapply(1,2)将参数1放到指定的方法2中
uni_glm <- lapply(variable.names,uni_glm_model)
uni_glm

#生成单变量分析的综合结果
library(plyr)
uni_glm <- ldply(uni_glm,data.frame)
uni_glm

#将单因素分析的 结果写道csv中
write.csv(uni_glm,"uni.csv")

#将p<0.05的结果挑选出来
uni_glm1 <- uni_glm[uni_glm$P<=0.05,]
uni_glm1





############################## 多因素分析 ######################################
#先写一个公式,方便后面重复利用
#collapse 用什么符号连接
fml <- as.formula(paste0('hypoglycemia==1~',paste0(uni_glm$characteristics[uni_glm$P<0.1],collapse = '+')))
fml

modelA <- glm(fml,data=dev,family = binomial)

modelX<-glm(hypoglycemia~1,data = dev,family=binomial)
modelX
#向前法取值
#step()逐步 scope范围 upper最多有几个,lower最少 family=binomial按向逻辑回归
#direction ="forward"方向向前
modelB<-step(modelX,scope=list(upper=~ course_of_disease + Hyperlipidemia + 
                                 Treat_Time + Education + gender + BUN + RBG + TC,
                               lower=~1),data = dev,family=binomial,direction ="forward")

summary(modelB)

#向后法取值backward
modelC <- step(modelA,direction ="backward")
summary(modelC)

#both法取值both
modelD<-step(modelA,direction = "both")
summary(modelD)



#看模型的系数及95%CI
cbind(coef=coef(modelD),confint(modelD))
#看模型的OR及95%CI
exp(cbind(OR=coef(modelD),confint(modelD)))
#模型AIC
AIC(modelA,modelB,modelC,modelD)

#模型AIC进行比较
anova(modelA,modelB,test = "Chisq")
anova(modelA,modelC,test = "Chisq")
anova(modelA,modelD,test = "Chisq")
anova(modelB,modelC,test = "Chisq")
anova(modelB,modelD,test = "Chisq")
anova(modelC,modelD,test = "Chisq")

#采用上述比较,决定用modelD的结果
modelD<-step(modelA,direction = "both")
modelD
glm3<-summary(modelD)
glm3

#将多因素分析结果制备为发表格式
glm3$coefficients

OR<-round(exp(glm3$coefficients[,1]),2)
OR
SE<-round(glm3$coefficients[,2],3)
CI2.5<-round(exp(coef(modelD)-1.96*SE),2)
CI97.5<-round(exp(coef(modelD)+1.96*SE),2)
CI<-paste0(CI2.5,'-',CI97.5)
B<-round(glm3$coefficients[,1],3)
Z<-round(glm3$coefficients[,3],3)
P<-round(glm3$coefficients[,4],3)

#制作数据框
mlogit <- data.frame(
  'B'=B,
  'SE'=SE,
  'OR'=OR,
  'CI'=CI,
  'Z' =Z,
  'P'=P)[-1,]
mlogit

#展示一下数据库变量列表,找到多因素分析的属性位置
names(mydata)
#提取最终模型,多因素分析有意义变量
multinames <- as.character(colnames(mydata)[c(4,6,9,13,30,32,38)])
multinames
#创建出characteristics列
mlogit<-data.frame('characteristics'=multinames,mlogit)



#合并先单后多分析表格
final <- merge.data.frame(uni_glm,mlogit,by = 'characteristics',all = T,sort = T)
#导出文件 
write.csv(final,"final.csv")





################################ 模型验证 ######################################
#上面获取到的模型
fml8<-as.formula(hypoglycemia == 1 ~ course_of_disease + Hyperlipidemia + Treat_Time + 
                   Education + gender + BUN + RBG + TC)
fml7<-as.formula(hypoglycemia == 1 ~ course_of_disease + Hyperlipidemia + Treat_Time + 
                   Education + BUN + RBG + TC)

model8<-glm(fml8,data = dev,family = binomial(logit))
model7<-glm(fml7,data = dev,family = binomial(logit))

#在建模人群中计算预测值
#predict()专门做预测值的函数 newdata对哪个数据集进行预测,用哪个模型,预测的值
dev$predmodel8<- predict(newdata=dev,model8,"response")
dev$predmodel7<- predict(newdata=dev,model7,"response")
View(dev)
#在验证人群中计算预测值
vad$predmodel8<- predict(newdata=vad,model8,"response")
vad$predmodel7<- predict(newdata=vad,model7,"response")
View(vad)





############################# 绘制ROC曲线(区分度) ##############################
#install.packages("pROC")
library(pROC)
#建模集,model8的auc与roc分析 smooth是否平滑
devmodelA <- roc(hypoglycemia~predmodel8,data=dev,smooth=F)
devmodelA

#提取devmodelA的auc保留3位(suc:曲线下的面积)
round(auc(devmodelA),3)
round(ci(auc(devmodelA)),3)#CI区间


#画图
#print.auc是否展示auc print.thres是否展示阈值 main标题名
plot(devmodelA, print.auc=TRUE, print.thres=TRUE,main = "ROC CURVE", 
     col= "blue",print.thres.col="blue",identity.col="blue",
     identity.lty=1,identity.lwd=1)

#数据集model7的auc与roc分析
devmodelB <- roc(hypoglycemia~predmodel7, data=dev,smooth = F)
round(auc(devmodelB),3)
round(ci(auc(devmodelB)),3)
plot(devmodelB, print.auc=TRUE, print.thres=TRUE,main = "ROC CURVE", 
     col= "red",print.thres.col="red",identity.col="red",
     identity.lty=1,identity.lwd=1)





############################### 绘制多条ROC曲线 ################################
#建模人群多条ROC方法一
#plot.roc(y,x,图的名字,百分位数,几号颜色)
#lines.roc()再添加一条roc曲线
devroc1 <- plot.roc(dev$hypoglycemia,dev$predmodel8,main ="dev ROC",percent=TRUE,col="1")
devroc2 <- lines.roc(dev$hypoglycemia,dev$predmodel7,percent=TRUE,col="2")
#图例
#legend(位置(右下角),图例名 = c("devmodel8","devmodel7"),颜色=c("1","2"),宽度 = 2)
legend("bottomright",legend = c("devmodel8","devmodel7"),col = c("1","2"),lwd = 2)

#建模人群roc比较
roc.test(devroc1,devroc2)





################################ 绘制校准度曲线 ################################
#方法一
install.packages("calibrate")
library(calibrate)
library(MASS)
install.packages("rms")
library(rms)

#在建模人群中绘制Calibration plot
#val.prob(p值,y)
val.prob(dev$predmodel8,dev$hypoglycemia)
val.prob(dev$predmodel7,dev$hypoglycemia)

#在建模人群中进行Hosmer-Lemeshow test检验
source("HLtest.R") #一定要把HLtest.R放在工作目录中
hl.ext2(dev$predmodel8,dev$hypoglycemia)
hl.ext2(dev$predmodel7,dev$hypoglycemia)
#在验证人群中进行Hosmer-Lemeshow test检验
source("HLtest.R") 
hl.ext2(vad$predmodel8,vad$hypoglycemia)
hl.ext2(vad$predmodel7,vad$hypoglycemia)



#方法二
install.packages("riskRegression")
library(riskRegression)
#7因子模型
formula<-hypoglycemia == 1 ~course_of_disease + Hyperlipidemia + Treat_Time + 
  Education + BUN + RBG + TC

#在建模集中制作校准曲线
fit1 <- glm(formula,data = dev,family = binomial )
#预测模型评分
xb <- Score(list("fit"=fit1),formula=hypoglycemia~1,#模型,y
            null.model = FALSE,                     #是否展现模型
            plots = c("calibration","ROC"),         #绘图需要什么
            metrics = c("auc","brier"),             #需要什么值
            B=1000,M=50,                            #内部验证
            data = dev)
#绘制校准曲线
plotCalibration(xb,col = "red")




################################# 绘制决策曲线 #################################
#重新读取数据
library(readr)
mydata <- read.csv("hypoglycemia_end.csv")
#删除有缺失值的行
mydata <- na.omit(mydata)
#数据集指定
dev = mydata[mydata$dataset==1,]
vad = mydata[mydata$dataset==0,]

#方法一(rmda验证完成)
#install.packages("rmda")
library(rmda)

#绘制DCA曲线
#8因子模型
model_1<-decision_curve(hypoglycemia ~course_of_disease + Hyperlipidemia + Treat_Time + 
                          Education + BUN + gender+ RBG + TC,
                        data = dev,
                        family = binomial(logit),#使用logit回归
                        thresholds = seq(0,1,by=0.01),#阈值0-1,以0.01为间隔
                        confidence.intervals = 0.95,#95%肯定区间ci值
                        study.design = 'case-control',#研究类型case-control
                        population.prevalence = 0.3)#粗使患病率300/1000
#7因子模型
model_2<-decision_curve(hypoglycemia ~course_of_disease + Hyperlipidemia + Treat_Time + 
                          Education + BUN + RBG + TC,
                        data = dev,
                        family = binomial(logit),
                        thresholds = seq(0,1,by=0.01),
                        confidence.intervals = 0.95,
                        study.design = 'case-control',
                        population.prevalence =0.3)

#绘制曲线
plot_decision_curve(model_1,curve.names = c('model_1'),
                    xlim = c(0,0.8),
                    cost.benefit.axis = FALSE,
                    col = c('red'),
                    confidence.intervals = FALSE,
                    standardize = FALSE)


plot_decision_curve(model_2,curve.names = c('model_2'),
                    xlim = c(0,0.8),
                    cost.benefit.axis = FALSE,
                    col = c('green'),
                    confidence.intervals = FALSE,    #TRUE
                    standardize = FALSE)

#绘制多条曲线
plot_decision_curve( list(model_1, model_2), 
                     curve.names = c("model_1", "model_2"),
                     col = c("blue", "red"), 
                     confidence.intervals = FALSE,  #remove confidence intervals
                     cost.benefit.axis = FALSE, #remove cost benefit axis
                     legend.position = "topright") #remove the legend "bottomright""none""topright"
#添加可信区间
plot_decision_curve( list(model_1, model_2), 
                     curve.names = c("model_1", "model_2"),
                     col = c("blue", "red"), 
                     confidence.intervals = TRUE,  #confidence intervals
                     cost.benefit.axis = FALSE, #remove cost benefit axis
                     legend.position = "topright") #add the legend "bottomright" "topright" "none"

#验证集决策曲线,需要先生成验证集的预测概率
fml8<-as.formula(hypoglycemia == 1 ~ course_of_disease + Hyperlipidemia + Treat_Time + 
                   Education + gender + BUN + RBG + TC)
fml7<-as.formula(hypoglycemia == 1 ~ course_of_disease + Hyperlipidemia + Treat_Time + 
                   Education + BUN + RBG + TC)


model8<-glm(fml8,data = dev,family = binomial(logit))
model7<-glm(fml7,data = dev,family = binomial(logit))

#在建模人群中计算预测值
dev$predmodel8<- predict(newdata=dev,model8,"response")
dev$predmodel7<- predict(newdata=dev,model7,"response")
#在验证人群中计算预测值
vad$predmodel8<- predict(newdata=vad,model8,"response")
vad$predmodel7<- predict(newdata=vad,model7,"response")

#8因子模型
vadmodel8 <- decision_curve(hypoglycemia~predmodel8,
                            data = vad,
                            fitted.risk = TRUE, 
                            thresholds = seq(0, .9, by = .05),
                            bootstraps = 200) #内部验证200次

plot_decision_curve(vadmodel8,curve.names = c('model_1'),
                    legend.position = "topright",
                    confidence.intervals = FALSE,    #remove confidence intervals)
                    standardize = FALSE) 


plot_decision_curve(vadmodel8,curve.names = c('model_1'),
                    legend.position = "topright",
                    confidence.intervals = TRUE,    #remove confidence intervals)
                    standardize = FALSE) 

#7因子模型
vadmodel7 <- decision_curve(hypoglycemia~predmodel7,
                            data = vad,
                            fitted.risk = TRUE, 
                            thresholds = seq(0, .9, by = .05),
                            bootstraps = 200) #内部验证200次

plot_decision_curve(vadmodel7, legend.position = "topright",
                    confidence.intervals = FALSE,
                    standardize = FALSE)  #remove confidence intervals

#绘制多条曲线
plot_decision_curve( list(vadmodel8, vadmodel7), 
                     curve.names = c("model_1", "model_2"),
                     col = c("blue", "red"), 
                     confidence.intervals = FALSE,  #remove confidence intervals
                     cost.benefit.axis = FALSE, #remove cost benefit axis
                     legend.position = "topright") #remove the legend "bottomright" "topright" "none"





#方法二(dca.R)
mydata <- read.csv("hypoglycemia_end.csv")
#删除有缺失值的行
mydata <- na.omit(mydata)
#数据集指定
dev = mydata[mydata$dataset==1,]
vad = mydata[mydata$dataset==0,]

source("dca.R")
#install.packages("nricens")
library(rms)
library(foreign)
library(nricens)

#创建2个回归模型用于演示
##8因素模型
modelA<-glm(hypoglycemia ~course_of_disease + Hyperlipidemia + Treat_Time + 
              Education + BUN + gender+ RBG + TC, data = dev, family = binomial(link="logit"),x=TRUE)
summary(modelA)
#用modelA预测概率(建模集和验证集)
dev$predmodelA<- predict(newdata=dev,modelA,"response")
vad$predmodelA<- predict(newdata=vad,modelA,"response")

##7因素模型
modelD <- glm(hypoglycemia ~course_of_disease + Hyperlipidemia + Treat_Time + 
                Education + BUN + RBG + TC, data = dev, family = binomial(link="logit"),x=TRUE)
summary(modelD)
#用modelD预测概率(建模集和验证集)
dev$predmodelD<- predict(newdata=dev,modelD,"response")
vad$predmodelD<- predict(newdata=vad,modelD,"response")

#训练集dca
dev <- as.data.frame(dev)#数据集转为数据框
dca(data=dev, outcome="hypoglycemia",
    predictors=c("predmodelA","predmodelD"),
    smooth="TRUE", probability=c("TRUE","TRUE"),
    xstop=0.5) 

#验证集dca
vad<-as.data.frame(vad)
dca(data=vad, outcome="hypoglycemia", 
    predictors=c("predmodelA","predmodelD"),
    smooth="TRUE", probability=c("TRUE","TRUE"),
    xstop=0.5) 

#







































































































